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Abstract. We discuss the current implementation of the ALI method 
into our HYDro dynamical RAdiation (Hydra) code for rapidly expanding, 
low density envelopes commonly found in core collapse and thermonuclear 
supernovae, novae and WR stars. Due to the low densities, non-thermal 
excitation by high energy photons (e.g. from radioactive decays) and 
the time dependence of the problem, significant departures from LTE are 
common throughout the envelope even at large optical depths. 

ALI is instrumental for both the coupling of the statistical equations 
and the hydrodynamical equations with the radiation transport (RT). 
We employ several concepts to improve the stability, and convergence 
rate/ control including the concept of leading elements, the use of net 
rates, level locking, reconstruction of global photon redistribution func- 
tions, equivalent-2-level approach, and predictive corrector methods. For 
appropriate conditions, the solution of the time-dependent rate equations 
can be reduced to the time-independent problem plus the (analytic) solu- 
tion of an ODE For the 3-D problem, we solve the radiation transport via 
the moment equations. To construct the Eddington tensor elements, we 
use a Monte Carlo scheme to determine the deviation of the solution for 
the RT equation from the diffusion approximation (ALI of second kind) . 

At the example of a subluminous, thermonuclear supernova (SN99by), 
we show an analysis of the light curves, flux and polarization spectra and 
discuss the limitations of our approach. 



1. Physical and Numerical Environment 



For supernovae, novae and Wolf Rayet stars, detailed hydrodynamical radiation 
calculations are required to provide a link between the observables such as light 
curves and spectra, and the underlying physics of the objects. These applica- 
tions go well beyond classical atmospheres. Density structures require detailed 
hydrodynamics, the low densities cause strong NLTE effects throughout the en- 
velopes, chemical profiles are depth dependent, the energy source and sink terms 
due to hydrodynamical effects and radioactive decays may dominate throughout 
the photon decoupling region, i.e. location of the photosphere and the physical 
properties are time dependent (see Fig. ||). Typically, velocity fields are of the 
order of 500 to 30,000 km/sec. Thus, as a major simplification, we can neglect 
the intrinsic line widths due to pressure broadening, magnetic fields, etc. 
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Figure 1. 

decay Ej, 



distance [10 1 / 6 cm] 

Temperature T, energy deposition due to radioactive 
Rosseland optical depth Tau (left scale) and density log(p) 
(right scale) are given as a function of distance (in 10 15 cm) for a typical 
SNe la at 15 days after the explosion. For comparison, we give the 
temperature T grey for the grey extended atmosphere. The light curves 
and spectra of Type la Supernovae are powered by energy release due 
to radioactive decay of 5e Ni — > 56 Co — > 56 Fe. The two dotted, vertical 
lines indicate the region of spectra formation. Most of the energy is 
deposited within the photosphere and, due to the small optical depth 
and densities, strong NLTE effects occur up to the very central region 
(from [Hoflich, 1995|). 



1.1. Numerical Tools and Methods 



The computational tools summarized below were used to carry out many of the 
analyzes of SNIa and Core Collapse Supernovae (Hoflich 1988, Hoflich, Miiller 
& Khokhlov 1993, [Hoflich, 1995|, [Howell et al. 2002|, ...)~^br the nomencla- 



ture we follow Mihalas & Mihalas (1988) if quantities are not defined explicitly. 
All components of the codes have been written or adopted in a modular form 
with well defined interfaces which allows an easy coupling (see Fig. |2| and code 
verification by exchanging e.g. various radiation transport modules but keeping 
identical the remaining setup (e.g. Fig. ||). The modules consist of physical 
units to provide a solution for e.g. the nuclear network, the statistical equations 
to determine the atomic level population, equation of states, the opacities, the 
hydro or the radiation transport problem. The individual modules are coupled 
explicitly. Consistency between the solutions is achieved iteratively by pertur- 
bation methods. Currently, not all modules can be combined simultaneously 
because not all iteration schemes have been implemented and because of re- 
quirements on CPU time: a) For full NLTE-spectra with large model atoms and 
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Numerical Environment for the ALI (Hydra-modules) 

Hydrodynamics (PPM) 



/ 



a) 1-D Lagrangian (spherical + front tracking) 

b) 3-D Eulerian (cartesian, PPM) 

c) Free expansion 



Nuclear network 

a) NSE 

b) Full network & decays 



EOS 

a) 1E10 > > 1 g/ccm 

b) lg/ccm > p 




MC gamma-ray transport 

a) 1-D spherical 

b) 3-D (given cartesian grid) 



LTE 



Statistical equations 
for ionization and 
level population 



Opacities 




Radiation transport (3 modules) 

al)Spherical, comoving Rybicki scheme (MKH75, 

76,81) for spherical LCs and atmospheres 
a2) Formal integration of RT in observes 
frame (spherical) 

b) Variable Eddington Tensor solver 
(implicit) for given Tensors 

bl) 1-D spherical (comoving) + energy 
b2) 3-D cartesian (observer) 

c) Monte Carlo Scheme 

cl) for Eddington tensor: 3-D,solve for difference 
between diffusion and R.T. equation ( ALI2) 
c2) Polarization: stationary transport 



Rem.: Not all modules can be combined simultaneously 

(Perturbation strategies and CPU-time: e.g. 3D-stmc.+ NLTE) 



Figure 2. Block diagram of our numerical scheme to solve radiation 
hydrodynamical problems including detailed equation of state, nuclear 
and atomic networks. For specific problems, a subset of the modules 
is employed (see text, and Figs. |I| & ||). The code is still under active 
development. For its current status, see the applications below. 



a high frequency resolution we are restricted to the time-independent case based 
on a given hydrodynamical structure (see Fig. |) and, for 3-D models, reduced 
atomic atomic models (level- merging/super- levels) have to be used (Fig. [|. b) 
Radiation hydrodynamics is restricted to a reduced frequency resolutions with 
reduced atomic levels (level-merging) and spherical geometry (see Fig. ||. In case 
of multi-dimensional radiation hydrodynamics, CPU-time requirements restrict 
applications even further. Currently, we can use of few frequencies to represent 
the fluxes e.g. in the Lyman and 'Balmer and higher continua', and 3-level 
atoms plus spherically symmetric velocity fields for the RT( [Hoflich, Khokhlov 
fc Wang 2002] ), or to the grey case for arbitrary field. Technically, we use the so 
called 'trivial' parallelization: Parallelization (via MPI and PVM) is achieved 
on the module basis, e.g. parallelization of the nuclear network over the grid 
points, using groups of photons for the Monte-Carlo radiation transport, and 
sub-domains with 'ghost-cells' with respect to the spacial and frequency coor- 
dinates for the hydro and the comoving frame radiation transport, respectively. 
In the remainder of this section, we want to describe the various modules. 
Hydrodynamics: The structure of the expanding envelopes are obtained by 
three different modules which by a) assuming free expansion, by solving the 
hydro equations in b) the Lagrangian frame for sp herical geometry including a 



front tracking scheme to resolve shock fronts (e.g. Hoflich, Wheeler fc Thiele- 
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mann 1998), or the Eulerian scheme for full 3-D using Cartesian coordinates 
based on PROMETHEUS ([Fryxefl, Arnett & Midler 199 1|). The hydro modules 



use an explicit Piecewise Parabolic Method (PPM) by |Collella &: Woodward 



(1984) to solve the compressible reactive flow equations. PPM is implemented 



as a step followed by separate remaps of the thermal and kinetic energy to avoid 
numerical generation of spurious pressure disturbances during propagation of 
reaction fronts (flames and detonations). 

Nuclear network &; high density EOS: Nuclear burning is taken into ac- 
count using Thielemann's nuclear reaction-network library ( [Thielemann, Nomoto 



& Hashimoto 1994, and referen ces therein). The main so urces for experimental 
rates are Caughlan et al.(1985) , Caughlan fc Fowler(1988) and Wiescher, Gorres] 



fc Thielemann (1990) . Typically, between 20 and 618 isotopes are taken into 
account. For the equation of state, we use a relativistic Fermi-gas with Coulomb 
corrections and crystallization, and radiation. 

Opacities and low density equation of state: For low densities, we solve the 
full rate equations (see below) to determine the level population, or assume LTE. 
Typically, about 500 to 600 discrete NLTE levels are included in full NLTE with 
about 20,000 to 40,000 line transitions. For the 3-D hydro, the effective number 
of levels is being reduced by about a factor of 10 by level-merging flHoflich 19901) , 
often referred to as the use of super-levels. Complete redistribution over each 
individual line both in frequency and in angle are assumed. This means that 
the relative populations within the sublevels or the merged levels are described 
by a Maxwell-Boltzmann distribution. The data for the atomic line transitions 



are taken from the new compilation of Kurucz (Kurucz 1991 , Kurucz 1995 ) and 
the opacity project (TOPBASE, e.g. |Cunto fe Mendoza 1992 ) from which about 



500,000 -2,000,000 lines are extracted, depending on the temperature and density 
range. Photon scattering by free electrons is included in the Klein-Nishina limit. 
Free-free cross sections are treated in the hydrogen approximation with free-free 
Gaunt factors according to |Seaton(1960)| and Gayet(1970)| . Radiative bound-free 
cross sections are taken from the opacity project (Z<28). 

Gamma-ray transport and energy deposition: The 7-ray transport is com- 
puted in spherical and th ree-dimension s using a Monte Carlo method (Hoflich, 
Khokhlov & Miiller 1992, Hoflich 2002 ) including relativistic effects and a con- 
sistent treatment of both the continua and line opacities. The interaction pro- 
cesses allowed are: Compton scattering according to the full angle-dependent 
Klein-Nishina formula, pair-production, and bound- free transitions (Ambwani 
fc Sutherland 1988| , prown fc Leventhal 1981 ). 

Radiation transport: Several modules are available. For the spherical geom- 
etry and in the stationary case, the radiation transport equation is solved in the 
comoving frame ( Mihalas, Kunacz &: Hummer 1975 , Mihalas, Kunacz fc Hum"] 
mer 1976| , Mihalas, Kunacz fc Hummer 1976b| , hereafter MKH- methods). For 



the time dependent cases, we use variable Eddington Tensor solvers (implicit 
in time) ( |Mihalas & Mihalas 19841 , [Stone, Mihalas & Norman 1992j , [Hoflich et| 
al. 1993] ). We assume that the Eddington factors are constant during each time 
step. In the spherical case, we use the MKH methods. To obtain the correct 
solution for the Eddington tensor elements in 3-D, we use a Monte Carlo method 
to compute the difference between the solution of the non-equilibrium diffusion 
and full radiation transport equation. We calculate the difference between the 
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solutions for computational accuracy and efficiency. In particular, the random 
element in MC avoids spatial fluctuations (and, thus, wiggles) in propagating, 
plane light fronts though, still, its speed in the free stream limit may differ from 
the light speed by 5 to 20 %. The same Monte Carlo solver is used which has 
been applied to compute 7-ray and polarization spectra for scattering dominated 
atmospheres (e.g. [Hoflich et al. 1992; , [Howell et al. 2002| ). The Monte Carlo 
method is appropriate for this problem because of its flexibility with respect to 
the geometrical and velocity structures. 

AutomaticMeshRefinement for radiation transport by a MC torch: 

Automatic Mesh Refinement is a well established procedure in hydrodynamics 
and allows to adjust the resolution to the requirements. The same holds true 
also for radiation transport. In stellar atmospheres, a logarithmic spacing of the 
optical depth r is adopted to guarantee an appropriate resolution. In dynami- 
cal problems, extended atmospheres/envelopes with arbitrary morphologies the 
problem boils down to determine the region of the photosphere, or the 'skin' of 
an optically thick object. We start with a grid, equally spaced in mass (in co- 
moving frame) or space (in Eulerian frame), and employ a 'Monte Carlo Torch' 
for the grid refinement. As a basic concept, the photosphere is defined as the 
region from which the photons can escape or, because of the symmetry of the 
problem, a photon can penetrate to when coming the outside. With the 'torch', 
we eluminate the object from the outside and calculate the path of the photons 
(with various energies) in random directions. If a photon is interacting in a 
computational cell, we increase its photon counter. Typically, the number of 
test photons is of the order of the number of cells npt in a grid times the number 
of representative frequencies u rep . In our examples (see Figs. || &i |5|), about 10 6 
to 10 7 are used both in the spherical (n r — 912 & Vrep — 

1000) ) and in the 3-D 

(n r sa 5 x 10 5 & 10 x v rep ). Computational cells are divided in half if their photon 
count exceeds the average by a factor of 5 to 10. By eluminating the object from 
the edges of the domain, 3-D surfaces will be traced only if they can be seen from 
the domain boundaries. Concave structures such as the inner edges of Rayleigh 
Taylor fingers will be missed. Therefore, ~ 10 additional photons are emitted 
from each grid cell. A photon triggers a counter only if it has crossed 2 or more 
cell boundaries to avoid rezoning of optically thick layers. The requirements of 
this approach depends on the physical situation. In our examples, we allow for 
rezoning every 20 to 50 time steps. To test the accuracy, we increase the number 
of test photons by a factor of 10 about every 10 rezoning steps. 



1.2. Accelerated Lambda Iteration & Friends 

The radiation field couples the local, statistical equations and the hydro equa- 
tions, resulting in complex systems of integro-differential equations of rank 
DIM = n{x)n(y)n{z)n(u) J2 e i Eion(eZ) T,ievd(el,ion) n (el, ion, level) where n(x,y,z) 
and n(i^) are the spatial and frequency coordinates, and n(el, ion, level) are the 
atomic level population of "level" for element "el" with the ionization "ion". 
The rank can easily exceed 10 12 -- 13 which prohibits a direct solution. How- 
ever, the iterative solution provides an effective way to separate global and local 
quantities. 

a) The method of accelerated lambda iteration is used to remove the global 
dependencies produced by the the radiation field. 
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b) The statistical equations and the energy terms in the hydro equations are 



solved by a partial linearization method (Mihalas & Mihalas 1984). When eval- 
uating the mean intensity, J„, in the rate equations two extreme cases are con- 
sidered in which the source function of the specific transition does or does not 
dominate the total source function. This approach further assumes that the 
variation of the total source function can be represented by the transition under 
consideration or by the change of the source function during the previous itera- 
tion step, respectively. The first case is quite similar to the assumptions of the 
equivalent two level approximation, but with the difference that the perturbation 
terms are used instead of the total rates. 

c) The concept of "leading" elements is introduced. An appropriate ordering 
of elements or groups of elements allows a separation of the equations, which 
are solved by partial linearization. With this technique, the system of equations 
remains small, independently of the problem. This property is needed to achieve 
numerical stability and computing efficiency. 

d) We use the concept of level-locking ( [Hoflich 1990| ), often called super-levels 
for numerical efficiency. For stability, we use level locking when the system is 
far from convergence. 

e) The scattering and thermal contributions to the source functions are sepa- 
rated by an equivalent two-level approximation for transitions from the ground 
levels. In effect, this allows to include the non-local nature of optically thick, 
scattering dominated transitions already during solution of the full radiation 
transport problem. In effect, this approach introduces an acceleration term for 



the convergence of ALI (Athay 1992, Avrett k Loeser 1992, Hoflich, 1995). 



f) Based on the explicit forms for the total opacity x{ v ) an d the source function 
S(u), we reconstruct redistribution functions in v for the photon, and limit the 
relative change between model iterations to ~ 10%. 

The accelerated lambda iteration: The intensity I can be obtained by the 
solution of the radiation transport equation 

^L + o P i = x(s-i) (i) 

where Op = Op(x, y, z,6/Sx, 5/5z, S/dz, 5/5v, and I is the intensity. On the right 
hand side, the source function S and and opacity x are local quantities. 

Formally, the solution of the radiation transport or, better, the first mo- 
mentum / can be written as J(x, y, z, v) = AS(x, y, z, t, 5/5x,5/5y, 5/6z), and J 
can be substituted in the rate and energy equation. The matrix elements Xji 
describe the contribution of S at the grid point j on J at the grid point i. In 
general, A is not constructed explicitly. Instead, it can be approximated by the 
following expression (e.g. Cannon 1973, S charmer, 1984| , Hillier 1990| , [Hoflich 



19901 , [Hubeny & Lanz 1992| , [Werner 1991] ) 



j(m) = A (m) s (m) _ A (m-1) s (m-l) + A * (m-1) ^(m) _ s (m-l) j 

and the solution is obtained iteratively where m stands for the iteration step. 
The first term corresponds to the 'classical' A iteration, which works well for 
low optical depths but the convergence of the A iteration becomes poor when A 
becomes diagonal and if the scattering part is large. The basic idea of the ALI is 
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introduce a second, term which removes the local contribution to the radiation 
field and drives the convergence. The implementation of the ALI method is 
reduced to the construction of A* and the way how the expression for J (eq. 2) 
is used for the rate and energy equations. 
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Figure 3. In expanding envelopes, a phoion travels both in the spa- 
tial and frequency space from (r ,u ) to (rx,v\) (red arrow). Conse- 
quently, any line transition (blue line) e[^ , ^i] will effect the absorp- 
tion probability of the photon whereas, in the static case, only lines 
will influence the absorption probability if \vu ne — v D \ is smaller than 
the internal line width Ai^j ne . For large velocities, a photon can be 
absorbed by a given line in small region with almost constant physical 
properties. This allows for an operator splitting of Op (eq. 1). 

The diagonal operator for rapidly expanding envelopes: In general, 
the spatial and frequency coordinates are coupled by the radiation transport 
equation, and the requirements on the corresponding resolutions are linked. In 
a Rybicki-like scheme for solving the radiation transport equation (e.g. Mihalas 
et al. 1975), the frequency spatial resolutions must be better the internal lines 
width, e.g. Au < vthermaU an d the velocity difference between two neighboring 
cells must be Av < Av/v*c with c being the speed of light. For SN, the velocities 
are of the order of 10 4 km/ sec, and the line width are given thermal velocities 
v therm ~ lkm/ sec which requires a > 10,000 grid points in each direction. 
However, the extreme ratio between line width and Doppler shift allows for an 
operator splitting of Op (eq. 1) into a directional and frequency dependent part 
(see Fig. ||). This allows reduce the solution of the RT problem into the two- 
stream approximation and to attribute the frequency shift to the opacity which 
removes the link between the required spatial and v resolution. 
Expansion opacities: In the presence of large velocity fields, the influence 
of the lines on the opacities can be well represented in the narrow line limit 
( ISobolev 1957| , |Castor 1974 ) fo r the given (non-LT E) level populations. We use 
expressions and presumptions ( Hoflich et al. 19"93| ) very similar to those but but 
for the comoving frame. The Sobolev opacity of a certain line 1 between a lower 
and upper level i and j, respectively, is represented by the expression 
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Figure 4. Comparison of a explosion model with the SNe Ia 1999by. 
We show spectra at day 11 (upper left), the the B and V light curves 
(right plots), and the chemical structure (lower panel). The explosion 
and evolution of the spectra are calculates self-consistently with the 
only free parameters being the initial structure of the exploding White 
Dwarf, and a parameterized description of the nuclear burning front. 
The explosion is calculated by a spherical hydro-modules using 912 
depth points including a nuclear network with 218 isotopes. After 
about 10 seconds, only radioactive decays are taken into account but 
we solve the time-dependent, full NLTE, radiation hydro. The LCs 
are based on several thousand NLTE-spectra utilizing 912 depth, 2000 
vs and atoms with a total of 50 super-levels. For several moments of 
time, 1, detailed NLTE spectra have been constructed using 90 depth & 

500 NLTE levels (from |H5fhch et al. 2002] ). 



~ 30, 000 v points with 

However, the probability to be absorbed by a given line needs to be modified 
by the probability that a photon did interact before with another line or a 
continuum. Similar to Karp et al.(1977) , we obtain the following expression if 



N scattering lines are ordered with respect to the wavelength 



with 
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Note that the Sobolev optical depth is formulated for a given reference 
frame. However, we can transform the opacity from frame to the frame of a 
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Figure 5. Comparison between theoretical spectra for the SNe la 
1999by at about maximum light assuming spherical symmetry. The 
spectra are based on our spherical (blue) and full 3-D radiation trans- 
port scheme, using 90 depth points, 20,000 frequency points and 520 
NLTE-levels and 67/67/67 depth points, 2000 frequency points and 50 
NLTE-super-levels, respectively. Differences are up to 20 %. They can 
be understood due to the lower resolution in the 3-D calculatio ns. As 
an application, asphericity effects of SN99by (see Fig. see Howell 
et gj? 2002|). 



neighboring grid point by 

x (v/c) = x(0)(l. 



+ 0(v/c) 



x(o) 



(6) 



if we assume a linear dependence of \u in the frequency space and peace-wise 



linear velocity fields (Hbflich 1990). Here, r is the local radius of curvature of 
the velocity field for a direction D. For A*, we can neglect the order 0(v/c) in 
the expression above. 

The matrix elements of the A* operator: We construct A* in 2-stream 
approximation, i.e. Op — ► 5/5D and use intensity like variables u = J + + J~ for 



direction D with the grid d. As in static atmospheres ( Olson, Auer fc Buchler 



1986|), we obtain 

5I_ 
Id 



x(I — S) and with 5t = x$d =3- —u = u — S (7) 

or 



By using the flux and intensity like variables we get the final equation 
Ud+l ( 1 1 \ , "d-i 



Ar d ■ Ar d+ i /2 



Ud 



+ 



+ 



Ar d -Ar d+1/2 Ar d • r d _ 1/2y / AT d _ 1/2 Ar d 



Ud ~ Sd- 
(8) 
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with 



S d = = u *d e W \~{ T d ~ Td-i)],u*d+i = «Sexp [-(to+i - T d )) (9) 



Xd±l/2 ■= 7, \Xd±l + Xd] Ar l+l/2 : = Xd+1/2 \ z d ~ Zd+l] 



(10) 



where u* d is the "flux" like quantity. After some transformations one gets the 
algebraic expression 



^ T d-i/2 ■= Xd-i/2 \ z d-\ ~ z d ] , Ar d+1/2 := - Ar d+1/2 + Ar d _!/ 



(11) 



/ 1 - exp(-Ar d ) 1 - exp(-Ar rf _i) 



Ar rf Ar d+1/2 



Ar d • Ar d _ 1/2 



(12) 



At the inner boundary of each ray we take for the optically thick and thin cores 
= (1 — e~ Ard = dmax ), respectively. From the the known 



d=dn 



1 * 

1 U d=dn 



we get the diagonal element by 



A? 



dfid(pu d . 



(13) 



One remaining problem related to the property that the locality of the oper- 
ator depends on the resolution, i.e. the convergence rate becomes grid dependent 
and, for high resolutions, the system may become unstable. The information on 
corrections propagates only one grid point per model interaction m. To over- 
come these problems and if the solution is far from convergence, we refine the 
the 'grid/average over several zones' corresponding to the thermalization depth 
T therm of the transition with the highest opacity of a given element. The es- 
timate of Ttherm based on the equivalent-two-level approach. Successively, we 
increase the effective resolution to the full grid resolution. For our applications, 
this approach has been found to be more stable rather than using off-diagonal 



elements (Werner 1991) but, in light of discussions at this meeting, the latter 
approach will be revisited. 

The matrix elements of the statistical equations: The matrix elements 
of the statistical equation A x = b of each species are given by 



-(Rij + Cij) i<j 
£fci(§)(^ + C u ) + Ei= i+ i(Rii + c u ) j = i 



(14) 



Ri j and C{ j are the radiative and collisional rates relating the levels i and j . 

The radiative rates between a lower and upper level i and j, respectively, 
are given by 

^^-Jvdv + R^-rays) (15) 
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e kr dv. 



(16) 
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Here, where the coupling between the radiation field and the statistical 
equations occurs, in principle, the solution of the statistical equations of the 
current model iteration (m) has to be known in order to calculate the new value 

of the radiation field jj^ . In order to estimate the updated rates, we distinguish 
two cases: 

Case A: S u /xu 3> Sij/xij and B) S u /xu ~ Sij/Xij- Case A: Eqs. 15 and 16 can 
be written as a sum of the "zero" order 



= 4vr| ^l^S^dv (17) 
otij{v) 



Rji = 47T 



c 2 



e kT d vdv (18) 



hv 

and by correction terms of "first" order 

RlY r = 4tt[ ^|Ma*(sM - Si m -^)dv (19) 

R-orr = 4n f ^M A .( S (m) _ 5 (m-l)) e ^ d v ( 20 ) 

J J hv 

In order to estimate the source functions we use again the concept of the 
"leading" elements. If the total emissivity at the considered frequency is domi- 
nated by an element which is higher in the hierarchy, the new difference between 
the old and the new source function is substituted by the difference of the dom- 
inating transitions (marked by the index dom), i.e. 

(t» " S^) - {Sfl ~ ■ X dom /Xv (21) 

Otherwise the monochromatic source function is just an extrapolation in 
the second order of the iteration steps, i.e. 

ASi m) ~ ASl m ~V + _ _ ^—^ + ... (22) 

with 

AS^ -=sW -S^l (23) 



Case B: If the transition in consideration dominates the total source function, 
we use the approach S u s=s Sij(u). In order to get expressions for the "first" order 
corrections to the rates, essentially, using the same assumptions as in the case 
of the equivalent two level approach. The explicit forms of the source functions 
allow us to substitute the departure coefficient of the lower level in the rate 
equations (e.g. Pauldrach 1987). After some simple transformations, one ends 
up with the expression for the first order correction of the rates 

Rt° rr = and BST = 4n f ^M A *>) ^ ~ S Jt^l dl/ (24) 

3 3 J hv yV ; Sp\v) 
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The separation between the cases A and B is needed to overcome numerical 
instabilities during the model iteration. The rate equations and the hydrody- 
namics are iteratively coupled via the energy equation (e.g. Hoflich et al. 1993| ). 



For J, we use the same expressions and case separation as just outlined. 
1.3. The Time-dependent Rate Equations 

An approximate solution of the time dependent non-LTE problem is outlined 
flHoflich 1991b ). The statistical equations are reduced to a form of time inde- 



pendent equations which can be solved by one of the standard methods for the 
time independent case, where the time dependent quantities enter the statistical 
equations just in an additional explicit rate. These are calculated analytically 
solving a linear inhomogeneous differential equation of first order. 

We assume that the time dependence is caused by a limited number of 
transitions. In our implementation, we assume that the time dependence is 
caused by the bound-free transitions. This condition is valid as long as forbidden 
transitions are suppressed because allowed bound-bound transitions are faster 
than bound-free transitions by orders of magnitude. This limits the applicability 
e.g., for SNela, to times < 40 to 60 days after the explosion. Forbidden lines 
causes the discrepancies in the late LCs in Fig. |j. 

In general, the time dependent rate equations are given by the expression 

■ 

+ V(riiu) = ^(rijPji - rtjPy) + n k (R ki + C ki ) - ni(R ik + C ik ) (25) 

where the rate Py is the sum of radiative and collisional processes. 

The time independent solution hi of the statistical equations is defined by 
the equations 

= ^(hjPji - hiPij) + h k (R ki + C ki ) - hi{R ik + C ik ) (26) 

that can be solved by a standard method. 

For most astrophysical applications, the time scales of allowed bound-bound 
transitions are much shorter than those of the bound-free transitions, and we 
obtain the relation 

Hi hi , . 

=> — = — . (27) 
rij raj- 
Subtracting the time independent rate equations from the time dependent 
ones both of which are normalized to and hi, respectively, results in 

i.e. an inhomogeneous linear differential equation of first order 

^^ = -C ini + C 2 (29) 



with 



Ufa 

Ci := —(R k i + C ki ) and C 2 := n k (R ki + C ki ). (30) 

Ui 
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This inhomogeneous ODE is solved analytically as an initial boundary problem. 
If more than two ionization stage are important, we end up with a system of 
ODEs to be solved numerically, and n~k are determined iteratively by using 
charge conservation in a simple fix point iteration, or a Newton-Raphson method. 

Acknowledgments. This research is supported in part by NASA Grant 
LSTA-98-022. 
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